---
title: "pa1_analyses"
output: word_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r load libraries}
library(readr)
library(reshape2)
library(plyr)
library(ggplot2)
library(mediation)
library(matrixStats)
library(tidyverse)

```

```{r import}

#load cleaned csv from matlab
fulldata <- read_csv("~/Desktop/Studies/PunishmentAversion/Study1/Final_Analyses/baseline/pa1_baseline_cleandata_withQualComp.csv",col_names = TRUE)

```

```{r clean_data}
                   
#Remove rows with any missing data
cleandata<-fulldata[complete.cases(fulldata), ]

#Specify that the variables are factors
cleandata$condition<-factor(cleandata$condition,labels=c("action","outcome", "punishment"))

#separate out the scenarios
actiondata=cleandata[which(cleandata$condition=="action"),]
outcomedata=cleandata[which(cleandata$condition=="outcome"),]
punishmentdata=cleandata[which(cleandata$condition=="punishment"),]
```

```{r calculate averages}
#calculate the average judgment for each punishment for each condition (saved as vectors)
start_col = 4
end_col <- 58

actionaverage = colMeans(actiondata[start_col:end_col])

outcomeaverage = colMeans(outcomedata[start_col:end_col])

punishmentaverage = colMeans(punishmentdata[start_col:end_col])

#create averages dataframe (combine vectors into single data frame)
averages.data <- data.frame(actionaverage, outcomeaverage,punishmentaverage)

#add scenario number as a variable 
averages.data$scenario=rownames(averages.data)

```
#Get relationship between action and outcome scores
```{r calculate averages}
#get correlation between action and outcome in full dataset
mosaic::cor.test(actionaverage ~ outcomeaverage, data = averages.data)

#plot correlation between action and outcome in full dataset
#get the model
act_out_mod <- lm(averages.data$actionaverage ~ averages.data$outcomeaverage, data=averages.data)
summary(act_out_mod)

#plot the model
plot(averages.data$actionaverage ~ averages.data$outcomeaverage, col="lightblue", pch=19, cex=2, data=averages.data,main = "Action vs Outcome", 
     xlab = "Average action judgment", ylab = "Average outcome judgment")

#add regression line
abline(act_out_mod, col="red", lwd=3)

#add labels to points (rownames)
text(averages.data$actionaverage ~ averages.data$outcomeaverage, labels=rownames(averages.data), data=averages.data, cex=0.9, font=2)

```

#Select subset of data to make action and outcome orthogonal
```{r subset data}

#create subset of scnearios that exclude those from the lower left and upper right corners of the action ~ outcome scatter plot
subset_new <- subset(averages.data, scenario==22|scenario==8|scenario==16|scenario==9|scenario==3|scenario==39|scenario==40|scenario==44|scenario==21|scenario==14|scenario==7|scenario==27|scenario==5|scenario==43|scenario==33|scenario==25|scenario==45|scenario==23|scenario==53|scenario==47|scenario==11|scenario==17|scenario==18|scenario==24|scenario==20|scenario==30|scenario==31|scenario==55|scenario==49|scenario==54|scenario==48|scenario==28|scenario==46)

#get number of scenarios in new subset
nrow(subset_new)

subset_new$scenario <- factor(subset_new$scenario, labels = c("hallucinate", "drug_for_life", "punch", "finger", "knock_out_teeth", "paralysis", "shock_no_effect", "slap", "poison_ill", "whip", "straightjacket", "iron_ball", "waterboard", "shirt", "siren", "labor_camp", "starve", "coal_miner", "toxic_waste", "house_arrest_20", "prison_1", "prison_5", "prison_20", "isolation_1", "isolation_5", "isolation_20", "notify_employers", "restrain_1", "restrain_5", "restrain_20", "national_news", "insults", "publish_diary"))

#output subset_new as a csv for averaging with replication data
write.csv(subset_new, "~/Desktop/Studies/PunishmentAversion/Study2/study1_baseline.csv", row.names=FALSE)

#get correlation between action and outcome in full dataset
mosaic::cor.test(actionaverage ~ outcomeaverage, data = subset_new)

# plot action vs outcome with new subset
aa <- subset_new$actionaverage
oa <- subset_new$outcomeaverage
pa <- subset_new$punishmentaverage

#get model of action predicted by outcome in new subset
act_out_mod2 <- lm(aa ~ oa, data=subset_new)
summary(act_out_mod2)

# plot the model
plot(aa ~ oa, col="lightblue", pch=19, cex=2, data=subset_new,main = "Action vs Outcome", 
     xlab = "Average action judgment", ylab = "Average outcome judgment")
#add regression line
abline(act_out_mod2, col="red", lwd=3)

#add labels to points (rownames)
text(aa ~ oa, labels=rownames(subset_new), subset_new, cex=0.9, font=2)

```
#Begin baseline analyses
#Baseline Plots
```{r basic_plots}
#plot action vs. outcome
ggplot(subset_new,aes(y=actionaverage,x=outcomeaverage))+
  geom_point()+
  geom_smooth(method=lm)


#plot punishment predicted by action
ggplot(subset_new,aes(y=punishmentaverage,x=actionaverage))+
  geom_point()+
  geom_smooth(method=lm)

#plot punishment predicted by outcome
ggplot(subset_new,aes(y=punishmentaverage,x=outcomeaverage))+
  geom_point()+
  geom_smooth(method=lm)

```

```{r baseline analysis}

# writing the full model
mod1 <- lm(pa~aa+oa)
summary(mod1)

#to test the hypothesis that the aa beta is significantly larger than the oa beta
library(car)
linearHypothesis(mod1, c("aa = oa"))
```

###INGROUP/OUTGROUP ANALYSES###
```{r import}

#load cleaned csv from matlab
ingroup_outgroupdata <- read_csv("~/Desktop/Studies/PunishmentAversion/Study1/Final_Analyses/ingroup_outgroup/pa1_ingroup_outgroup_cleandata_withQualComp.csv", col_names = TRUE)

```
#Clean ingroup/outgroup data
```{r clean_data}

#Remove rows with any missing data
clean_inout_data <- ingroup_outgroupdata[complete.cases(ingroup_outgroupdata), ]

#Remove any subject who did not pass the qualitative comprehension
clean_inout_data <- clean_inout_data[which(clean_inout_data$qual_comp == 1), ]

#Specify that the variables are factors
clean_inout_data$condition<-factor(clean_inout_data$condition,labels=c("ingroup","outgroup"))

#separate out the scenarios
ingroup_punishmentdata=clean_inout_data[which(clean_inout_data$condition=="ingroup"),]
outgroup_punishmentdata=clean_inout_data[which(clean_inout_data$condition=="outgroup"),]
```

#Calculate ingroup/outgroup averages
```{r calculate averages}
#get the column numbers for the first and last scenario columns
first_col <- which(colnames(clean_inout_data) == "shirt")
last_col <- which(colnames(clean_inout_data) == "restrain_20")

#Calculate the average judgment for each punishment for each condition (saved as vectors)
ingroup_punishmentaverage = colMeans(ingroup_punishmentdata[first_col:last_col])
outgroup_punishmentaverage = colMeans(outgroup_punishmentdata[first_col:last_col])

#create final averages dataframe make 2 and join them together
ingroupdata <- data.frame(subset_new$actionaverage, subset_new$outcomeaverage,ingroup_punishmentaverage)
colnames(ingroupdata)[1] <- c("actionaverage")
colnames(ingroupdata)[2] <- c("outcomeaverage")
colnames(ingroupdata)[3] <- c("punishmentaverage")
#set the condition to 1 when it is ingroup condition
ingroupdata$condition = 1

#create a variable that specifies the scenario
ingroupdata <- rownames_to_column(ingroupdata, var = "scenario")

outgroupdata <-data.frame(subset_new$actionaverage, subset_new$outcomeaverage,outgroup_punishmentaverage)
colnames(outgroupdata)[1] <- c("actionaverage")
colnames(outgroupdata)[2] <- c("outcomeaverage")
colnames(outgroupdata)[3] <- c("punishmentaverage")

#set the condition to -1 when it is outgroup
outgroupdata$condition = -1

#create a variable that specifies the scenario
outgroupdata <- rownames_to_column(outgroupdata, var = "scenario")

#combine into a full dataset
inout_averages.data <- rbind(ingroupdata,outgroupdata)

#mean center numerical variables
inout_averages.data$actionavg_cent<-inout_averages.data$actionaverage-(mean(inout_averages.data$actionaverage))
inout_averages.data$outcomeavg_cent<-inout_averages.data$outcomeaverage-(mean(inout_averages.data$outcomeaverage))
inout_averages.data$punishmentavg_cent<-inout_averages.data$punishmentaverage-(mean(inout_averages.data$punishmentaverage))

#create interaction variables
inout_averages.data$aa_x_cond <- as.numeric(inout_averages.data$condition)*inout_averages.data$actionavg_cent
inout_averages.data$oa_x_cond <- as.numeric(inout_averages.data$condition)*inout_averages.data$outcomeavg_cent

#specify that categorical variables are factors
inout_averages.data$scenario <- factor(inout_averages.data$scenario)
inout_averages.data$condition <- factor(inout_averages.data$condition, labels = c("outgroup", "ingroup"))

write.csv(inout_averages.data, "~/Desktop/Studies/PunishmentAversion/Study2/study1_inout.csv", row.names=FALSE)
```
#Ingroup/outgroup descriptives
```{r descriptives}
punishment_desc <- ddply(inout_averages.data, c("condition"), summarise,
               N    = length(punishmentaverage),
               mean = mean(punishmentaverage),
               sd   = sd(punishmentaverage),
               se   = sd / sqrt(N)
)
print(punishment_desc)

ingroup_data<-inout_averages.data[which(inout_averages.data$condition=="ingroup"),]
outgroup_data<-inout_averages.data[which(inout_averages.data$condition=="outgroup"),]

mosaic::cor.test(punishmentavg_cent~actionavg_cent,data=ingroup_data)
mosaic::cor.test(punishmentavg_cent~actionavg_cent,data=outgroup_data)

mosaic::cor.test(punishmentavg_cent~outcomeavg_cent,data=ingroup_data)
mosaic::cor.test(punishmentavg_cent~outcomeavg_cent,data=outgroup_data)
```

#Ingroup/outgroup analyses
```{r linear models}

# full model
full_mod <- lme4::lmer(punishmentaverage ~ actionavg_cent + outcomeavg_cent + condition + aa_x_cond + oa_x_cond + (1|scenario),data=inout_averages.data)

summary(full_mod)

#to test the hypothesis that the actionavg beta is significantly larger than the outcomeavg beta
library(car)
linearHypothesis(full_mod, c("actionavg_cent = outcomeavg_cent"))

#model without action main effect
NoAction_mod <- lme4::lmer(punishmentaverage ~ outcomeavg_cent + condition + aa_x_cond + oa_x_cond + (1|scenario),data=inout_averages.data)

summary(NoAction_mod)
anova(full_mod,NoAction_mod)

#model without outcome main effect
NoOutcome_mod <- lme4::lmer(punishmentaverage ~ actionavg_cent + condition + aa_x_cond + oa_x_cond + (1|scenario),data=inout_averages.data)

summary(NoOutcome_mod)
anova(full_mod,NoOutcome_mod)

#model without condition main effect
NoCondition_mod <- lme4::lmer(punishmentaverage ~ actionavg_cent + outcomeavg_cent + aa_x_cond + oa_x_cond + (1|scenario),data=inout_averages.data)

summary(NoCondition_mod)
anova(full_mod,NoCondition_mod)

mosaic::favstats(punishmentaverage ~ condition, data = inout_averages.data)

#model without action x condition interaction
NoActxCond_mod <- lme4::lmer(punishmentaverage ~ actionavg_cent + outcomeavg_cent + condition + oa_x_cond + (1|scenario),data=inout_averages.data)

summary(NoActxCond_mod)
anova(full_mod,NoActxCond_mod)

#model without outcome x condition interaction
NoOutxCond_mod <- lme4::lmer(punishmentaverage ~ actionavg_cent + outcomeavg_cent + condition + aa_x_cond + (1|scenario),data=inout_averages.data)

summary(NoOutxCond_mod)
anova(full_mod,NoOutxCond_mod)

```
#Ingroup/outgroup plots
```{r new graphs}

# action vs punishment + outcome vs punishment divided by condition

# punishment vs action 
ggplot(inout_averages.data,aes(y=punishmentaverage,x=actionaverage))+
  geom_point()+
  facet_wrap(~condition)+
  geom_smooth(method=lm)

# punishment vs outcome 
ggplot(inout_averages.data,aes(y=punishmentaverage,x=outcomeaverage))+
  geom_point()+
  facet_wrap(~condition)+
  geom_smooth(method=lm)

#action vs. outcome
ggplot(inout_averages.data,aes(y=actionaverage,x=outcomeaverage))+
  geom_point()+
  geom_smooth(method=lm)

#bar graphs

#punishment
ggplot(data=inout_averages.data, aes(x=scenario, y=punishmentaverage, fill=condition)) +
  geom_bar(stat="identity", width=.9, position=position_dodge()) + scale_fill_brewer(palette="Set2")

```

###SELF/OTHER ANALYSES###

```{r import}

#load cleaned csv from matlab
self_otherdata <- read_csv("~/Desktop/Studies/PunishmentAversion/Study1/Final_Analyses/self_other/pa1_self_other_cleandata_withQualComp.csv", col_names = TRUE)

```
#Clean self/other data
```{r clean_data}

#Remove rows with any missing data
clean_selfother_data <- self_otherdata[complete.cases(self_otherdata), ]

#Remove any subject who did not pass the qualitative comprehension
clean_selfother_data <- clean_selfother_data[which(clean_selfother_data$qual_comp == 1), ]

#Specify that the variables are factors
clean_selfother_data$condition<-factor(clean_selfother_data$condition,labels=c("self","other"))

#separate out the scenarios
self_punishmentdata=clean_selfother_data[which(clean_selfother_data$condition=="self"),]
other_punishmentdata=clean_selfother_data[which(clean_selfother_data$condition=="other"),]
```

#Calculate self/other averages
```{r calculate averages}
#get the column numbers for the first and last scenario columns
first_col <- which(colnames(clean_selfother_data) == "shirt")
last_col <- which(colnames(clean_selfother_data) == "restrain_20")

#Calculate the average judgment for each punishment for each condition (saved as vectors)
self_punishmentaverage = colMeans(self_punishmentdata[first_col:last_col])
other_punishmentaverage = colMeans(other_punishmentdata[first_col:last_col])

#create final averages dataframe make 2 and join them together
selfdata <- data.frame(subset_new$actionaverage, subset_new$outcomeaverage,self_punishmentaverage)
colnames(selfdata)[1] <- c("actionaverage")
colnames(selfdata)[2] <- c("outcomeaverage")
colnames(selfdata)[3] <- c("punishmentaverage")
#set the condition to 1 when it is self condition
selfdata$condition = 1

#create a variable that specifies the scenario
selfdata <- rownames_to_column(selfdata, var = "scenario")

#create final averages dataframe make 2 and join them together
otherdata <- data.frame(subset_new$actionaverage, subset_new$outcomeaverage,other_punishmentaverage)
colnames(otherdata)[1] <- c("actionaverage")
colnames(otherdata)[2] <- c("outcomeaverage")
colnames(otherdata)[3] <- c("punishmentaverage")
#set the condition to 1 when it is other condition
otherdata$condition = -1

#create a variable that specifies the scenario
otherdata <- rownames_to_column(otherdata, var = "scenario")

#combine into a full dataset
selfother_averages.data <- rbind(selfdata,otherdata)

#mean center numerical variables
selfother_averages.data$actionavg_cent<-selfother_averages.data$actionaverage-(mean(selfother_averages.data$actionaverage))
selfother_averages.data$outcomeavg_cent<-selfother_averages.data$outcomeaverage-(mean(selfother_averages.data$outcomeaverage))
selfother_averages.data$punishmentavg_cent<-selfother_averages.data$punishmentaverage-(mean(selfother_averages.data$punishmentaverage))

#create interaction variables
selfother_averages.data$aa_x_cond <- as.numeric(selfother_averages.data$condition)*selfother_averages.data$actionavg_cent
selfother_averages.data$oa_x_cond <- as.numeric(selfother_averages.data$condition)*selfother_averages.data$outcomeavg_cent

#specify that categorical variables are factors
selfother_averages.data$scenario <- factor(selfother_averages.data$scenario)
selfother_averages.data$condition <- factor(selfother_averages.data$condition, labels = c("other", "self"))

write.csv(selfother_averages.data, "~/Desktop/Studies/PunishmentAversion/Study2/study1_selfother.csv", row.names=FALSE)

```

#Self/Other descriptives
```{r descriptives}
punishment_desc <- ddply(selfother_averages.data, c("condition"), summarise,
               N    = length(punishmentaverage),
               mean = mean(punishmentaverage),
               sd   = sd(punishmentaverage),
               se   = sd / sqrt(N)
)
print(punishment_desc)

self_data<-selfother_averages.data[which(selfother_averages.data$condition=="self"),]
other_data<-selfother_averages.data[which(selfother_averages.data$condition=="other"),]

mosaic::cor(punishmentavg_cent~actionavg_cent,data=self_data)
mosaic::cor(punishmentavg_cent~actionavg_cent,data=other_data)

mosaic::cor(punishmentavg_cent~outcomeavg_cent,data=self_data)
mosaic::cor(punishmentavg_cent~outcomeavg_cent,data=other_data)
```

#Self/Other analyses
```{r linear models}

# full model
full_mod <- lme4::lmer(punishmentavg_cent ~ actionavg_cent + outcomeavg_cent + condition + aa_x_cond + oa_x_cond + (1|scenario),data=selfother_averages.data)

summary(full_mod)

#to test the hypothesis that the actionavg beta is significantly larger than the outcomeavg beta
library(car)
linearHypothesis(full_mod, c("actionavg_cent = outcomeavg_cent"))

#model without action main effect
NoAction_mod <- lme4::lmer(punishmentavg_cent ~ outcomeavg_cent + condition + aa_x_cond + oa_x_cond + (1|scenario),data=selfother_averages.data)

summary(NoAction_mod)
anova(full_mod,NoAction_mod)

#model without outcome main effect
NoOutcome_mod <- lme4::lmer(punishmentavg_cent ~ actionavg_cent + condition + aa_x_cond + oa_x_cond + (1|scenario),data=selfother_averages.data)

summary(NoOutcome_mod)
anova(full_mod,NoOutcome_mod)

#model without condition main effect
NoCondition_mod <- lme4::lmer(punishmentavg_cent ~ actionavg_cent + outcomeavg_cent + aa_x_cond + oa_x_cond + (1|scenario),data=selfother_averages.data)

summary(NoCondition_mod)
anova(full_mod,NoCondition_mod)

#model without action x condition interaction
NoActxCond_mod <- lme4::lmer(punishmentavg_cent ~ actionavg_cent + outcomeavg_cent + condition + oa_x_cond + (1|scenario),data=selfother_averages.data)

summary(NoActxCond_mod)
anova(full_mod,NoActxCond_mod)

#model without outcome x condition interaction
NoOutxCond_mod <- lme4::lmer(punishmentavg_cent ~ actionavg_cent + outcomeavg_cent + condition + aa_x_cond + (1|scenario),data=selfother_averages.data)

summary(NoOutxCond_mod)
anova(full_mod,NoOutxCond_mod)

```

#Self/Other plots
```{r new graphs}

# action vs punishment + outcome vs punishment divided by condition

# punishment vs action 
ggplot(selfother_averages.data,aes(y=punishmentaverage,x=actionaverage))+
  geom_point()+
  facet_wrap(~condition)+
  geom_smooth(method=lm)

# punishment vs outcome 
ggplot(selfother_averages.data,aes(y=punishmentaverage,x=outcomeaverage))+
  geom_point()+
  facet_wrap(~condition)+
  geom_smooth(method=lm)

#action vs. outcome
ggplot(selfother_averages.data,aes(y=actionaverage,x=outcomeaverage))+
  geom_point()+
  geom_smooth(method=lm)

#bar graphs

#punishment
ggplot(data=selfother_averages.data, aes(x=scenario, y=punishmentaverage, fill=condition)) +
  geom_bar(stat="identity", width=.9, position=position_dodge()) + scale_fill_brewer(palette="Set2")

#model without action or action x condition interaction
NoAct_NoActxCond.mod <- lme4::lmer(punishmentaverage ~ outcomeavg_cent + condition + oa_x_cond + (1|scenario),data=selfother_averages.data)
#get residuals from this model
selfother_averages.data$act_resids<-resid(NoAct_NoActxCond.mod)
#plot residuals predicted by action and condition
ggplot(selfother_averages.data,aes(y=act_resids,x=actionaverage))+
  geom_point()+
  facet_wrap(~condition)+
  geom_smooth(method=lm)

#model without outcome or outcome x condition interaction
NoOut_NoOutxCond.mod <- lme4::lmer(punishmentaverage ~ actionavg_cent + condition + aa_x_cond + (1|scenario),data=selfother_averages.data)
#get residuals from this model
selfother_averages.data$out_resids<-resid(NoOut_NoOutxCond.mod)
#plot residuals predicted by action and condition
ggplot(selfother_averages.data,aes(y=out_resids,x=outcomeaverage))+
  geom_point()+
  facet_wrap(~condition)+
  geom_smooth(method=lm)

```